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Abstract 

This Mathematica 7.0/8.0 package upgrades and extends the quantum computer 
simulation code called QDENSITY. Use of the density matrix was emphasized in 
QDENSITY, although that code was also applicable to a quantum state description. 
In the present version, the quantum state version is stressed and made amenable to 
future extensions to parallel computer simulations. The add-on QCWAVE extends 
QDENSITY in several ways. The first way is to describe the action of one, two and 
three- qubit quantum gates as a set of small (2 x 2, 4 x 4 or 8 x 8) matrices acting on 
the 2""? amplitudes for a system of Uq qubits. This procedure was described in our 
parallel computer simulation QCMPI and is reviewed here. The advantage is that 
smaller storage demands are made, without loss of speed, and that the procedure can 
take advantage of message passing interface (MPI) techniques, which will hopefully 
be generally available in future Mathematica versions. 

Another extension of QDENSITY provided here is a multiverse approach, as de- 
scribed in our QCMPI paper. This multiverse approach involves using the present 
slave-master parallel processing capabilities of Mathematica 7.0/8.0 to simulate er- 
rors and error correction. The basic idea is that parallel versions of QCWAVE run 
simultaneously with random errors introduced on some of the processors, with an 
ensemble average used to represent the real world situation. Within this approach, 
error correction steps can be simulated and their efficacy tested. This capability 
allows one to examine the detrimental effects of errors and the benefits of error 
correction on particular quantum algorithms. 

Other upgrades provided in this version includes circuit-diagram drawing com- 
mands, better Dirac form and amplitude display features. These are included in the 
add-ons QCWave.m and Circuits. m, and are illustrated in tutorial notebooks. 

In separate notebooks, QCWAVE is applied to sample algorithms in which the 
parallel multiverse setup is illustrated and error correction is simulated. These ex- 
tensions and upgrades will hopefully help in both instruction and in application to 
QC dynamics and error correction studies. 
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Program Summary 



Title of program: QCWAVE. Catalogue identifier: 



Program summary http://cpc.cs.qub.ac.uk/summaries 



Program available from: CPC Program Library, Queen's University of Belfast, 
N. Ireland. 

Operating systems: Any operating system that supports Mathematica; tested 
under Microsoft Windows XP, Macintosh OSX, and Linux FC4. 
Programming language used: Mathematica 7.0. 

Number of bytes in distributed program, including test code and documenta- 
tion: XX 

Distribution format: tar.gz 

Nature of Problem: Simulation of quantum circuits, quantum algorithms, noise 
and quantum error correction. 

Method of Solution: A Mathematica package containing commands to create 
and analyze quantum circuits is upgraded and extended, with emphasis on 
state amplitudes. Several Mathematica notebooks containing relevant exam- 
ples are explained in detail. The parallel computing feature of Mathematica is 
used to develop a multiverse approach for including noise and forming suitable 
ensemble averaged density matrix evolution. Error correction is simulated. 
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1 INTRODUCTION 



In this paper, QDENSITY [T] (a Mathematica |2] package that provides a 
flexible simulation of a quantum computer) is extended and upgraded by an 
add-on called QCWAVEQ The earlier flexibility in QDENSITY is enhanced 
by adopting a simple state vector approach to initializations, operators, gates, 
and measurements. Although the present version stresses a state vector ap- 
proach the density matrix can always be constructed and examined. Indeed, a 
parallel universe (or multiverse) approach is also included, using the present 
Mathematica 7.0/8.0 slave-master feature. This multiverse approach, which 
was published [5] in our QCMPI paper[^ allows separate dynamical evolu- 
tions on several processors with some evolutions subject to random errors. 
Then an ensemble average is performed over the various processors to pro- 
duce a density matrix that describes a QC system with realistic errors. Error 
correction methods can also be invoked on the set of processors to test the 
efficacy of such methods on selected QC algorithms. 



In section [2} we introduce qubit state vectors and associated amplitudes for 
one, two and multi-qubit states. In section |3j a method for handling one, 
two and three- qubit operators acting on state vectors with commands from 
QCWAVE are presented. 

In section |4| illustrations of how to apply gates to states are shown. In section 
[5| the multiverse approach is described and the parallel method for introduc- 
tion of errors and error correction are given. The ensemble averaged density 
matrix is then constructed. 



Additional upgrades, such as Dirac and amplitude displays and circuit draw- 
ing are presented in section |6| Suggested applications are presented in the 
conclusion section [7l 



^ Other authors have also developed Mathematica/Maxima QDENSITY based 
quantum computing simulations |3|4j . Hopefully, they will incorporate the ideas 
we provide herein in their future efforts. 

^ QCMPI is a quantum computer (QC) simulation package written in Fortran 90 with 
parallel processing capabilities. 
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2 MULTI-QUBIT STATES 



2.1 One-qubit states 

The basic idea of a quantum state, its representation in Hilbert space and the 
concepts of quantum computing have been discussed in many texts [61I71I8] . 
A brief review was given in our earher papers in this series [HIS]. Here we 
proceed directly from one, two and muhi-qubit states and their amphtudes to 
how various operators aher those amphtudes. 

To start, recaU that when one focuses on just two states of a quantum system, 
such as the spin part of a spin-1/2 particle, the two states are represented 
as either | 0) or | 1). A one qubit state is a superposition of the two states 
associated with the above and 1 bits: 

|vI/i) = Co|0)+Ci|l), (1) 

where Cq = (0 | \E'i) and Ci = (1 | \E'i) are complex probability amplitudes 
for finding the qubit in the state | 0) or | 1), respectively. The normalization 
of the state | \E'i) = 1, yields | Co P + I C*! P= 1- Note that the spatial 
aspects of the wave function are being suppressed; which corresponds to the 
particle being in a fixed location. The kets | 0) and | 1) can be represented 
as I 0) (q) and | 1) — )■ (5) • Hence a 2 x 1 matrix representation of this 
one-qubit state is: | ^i) ~^ (ci) • 

An essential point is that a quantum mechanical (QM) system can exist in 
a superposition of these two bits; hence, the state is called a quantum-bit 
or "qubit". Although our discussion uses the notation of a system with spin 
1/2, it should be noted that the same discussion applies to any two distinct 
quantum states that can be associated with | 0) and | 1). 

2.2 Two-qubit states 

The single qubit case can now be generalized to multiple qubits. Consider 
the product space of two qubits both in the "up" | 0) state and denote that 
product state as | 0) =| 0) | 0), which clearly generalizes to 

I Qi (I2) =\qi) I 92), (2) 

where gi,g2 in general take on the values and 1. This product is called a 
tensor product and is symbolized as 

I ?i 92) =1 qi)® I qi)- (3) 
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In QDENSITY, the kcts | 0), | 1) are invoked by the commands Ket[0] and 
Ket[l], and the product state by for example | 00) = Ket[0] Ket[0]. 

The kets | 00), | 01), | 10), & | 11) can be represented as 4 x 1 matrices 
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1 00) ^ 




;|01)^ 




; 1 10) ^ 




















1 























Hence, a 4 x 1 matrix 


representation of the two-qubit state 




1 * 


2) = 


Co 1 00) + Ci 


01) + C2 1 


10) + Cs 1 11), 





(4) 



(5) 



is: 



*2) ^ 



C2 



(6) 



Again Co = (00 | ^2),Ci = (01 | ^2),C2 = (10 | ^^2), and C3 = (11 | ^2), are 
complex probability amplitudes for finding the two-qubit system in the states 
I qi 52)- The normalization of the state (^^2 | ^2) = 1, yields 

I Co I' + I Ci |2 + I C2 p + I C3 1. (7) 



Note that wc label the amplitudes using the decimal equivalent of the bit 
product qi q2, so that for example a binary label on the amplitude Cio is 
equivalent to the decimal label C2. 



2.3 MuUi-qubit states 



For Uq qubits the computational basis of states generalizes to: 

I ^)n, =\qi)---\ (In,) =\(ll(l2 ■■■ Qn,) =| Q)- (8) 

We use the convention that the most significant qubit is labeled as qi and the 
least significant qubit by g^^. Note wc use q^ to indicate the quantum number 
of the ith qubit. The values assumed by any qubit is limited to either qi = 
or 1. The state label Q denotes the qubit array Q = (^i, ^2, • " " > Quq) , which 
is a binary number label for the state with equivalent decimal label n. This 
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decimal multi-qubit state label is related to the equivalent binary label by 

Uq 

n = qi- 2""-^ + q2 ■ 2"^-^ + . . . + . 2^ q, ■ T^"'' . (9) 

1=1 

Note that the ith qubit contributes a value of qi ■ 2"''"* to the decimal number 
n. Later we will consider "partner states" (| no), | ni)) associated with a given 
n, where a particular qubit is has a value of g^^ — 0, 

no = n - g,, • 2"''-^^ (10) 

or a value of qi^ — 1, 

= -l)-2"«"''- (11) 

These partner states are involved in the action of a single operator acting on 
qubit is, as described in the next section. 

A general state with Uq qubits can be expanded in terms of the above compu- 
tational basis states as follows 

l*)n,=E^QlQ)^ E Cn \n), (12) 

Q n=0 

where the sum over Q is really a product of Ug summations of the form J2qi=o,i ■ 
The above Hilbert space expression maps over to an array, or column vector, 
of length 2"« 



Co 



G 



L'0---00 



or with binary labels 



Co.. 



■01 



Ci., 



■11 



(13) 



The expansion coefficients C„ (or Cq) are complex numbers with the physical 
meaning that C„ = (n | ^)ng is the probability amplitude for finding the 
system in the computational basis state | n), which corresponds to having the 
qubits pointing in the directions specified by the binary array Q. Switching 
between decimal n and equivalent binary Q labels is accomplished by the 
Mathematica command IntegerDigits. 



In general, the complex amplitudes C„ vary with time and are changed by the 
action of operators or gates, as outlined next. 
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3 MULTI-QUBIT OPERATORS 



Operators that act in the multi-qubit space described above can be generated 
from a set of separate Pauh operators |3] acting in each qubit space. These 
separate Pauh operators refer to distinct quantum systems and hence they 
commute. Note, Pauh operators acting on the same qubit do not commute; 
indeed, they have the property a^aj — ajai = 2 i eijk cr^. The Pauh operator (Tq 
is just the unit 2x2 matrix. We denote a Pauh operator acting on qubit is as 
(tI^''\ where k = {x, y, z) = (1, 2, 3) is the component of the Pauh operator. For 
example, the tensor product of two qubit operators has the foUowing structure 



(ai I cxj I 6i)(a2 | aj \ 62) = (0-102 | (^i^^'^f'^ \ 

= (aia2 I ^ ® af \ hh) , (14) 

which defines what we mean by the tensor product of two qubit operators 
. The generahzation to more qubits is immediate 



3.1 One-qubit operators 



One-qubit operators change the amphtude coefficients of the quantum state. 
The NOT and Hadamard "H are examples of one-qubit operators of particu- 
lar interest: NOT = cr^ = ( ? J ) , 7/ = = -i^ ( } „\ ) . These have the 

following effect on the basis states NOT | 0) =| 1), NOT | 1) =| 0), and 
7^|0) = M),andH|l) = M). 

General one-qubit operators can also be constructed from the Pauli operators; 
we denote the general one-qubit operator acting on qubit s as fig- Consider 
the action of such an operator on the multi-qubit state | : 



^s|^k=ECQ ^^.IQ) (16) 
Q 

= J2 ■■■ ■■■ 111)--- {^s\qs)) - --\qn,)- 



91=0,1 q,=0,l g„„=0,l 



(17) 



^ The Pauli operators act in the qubit Hilbert space, and have the matrix repre- 
sentation: = ( ? J ) ; cTy = ( 5 ) ; cr^ = ( J ) . Here I = . 
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Here Qg is assumed to act only on the qubit is of value g^. The {Qs\ Qs)) term 
can be expressed as 



9^=0,1 



(18) 



using the closure property of the one qubit states. Thus Eq. (17) becomes 



^)n,=J2CQ^s\ Q) = 

Q 

■■■ ■■■ Y'^Q (9« I 1 9.) I ?i) ■ ■ ■ I gl) ■ ■ ■ kn,). 

91=0,1 gs=0,l g„„=0,l 9^=0,1 



(19) 



Now we can interchange the labels Qs ^ q's, and use the label Q to obtain the 
algebraic result for the action of a one-qubit operator on a multi-qubit state 



ns\^)n, = YCQ IQ)= E Cn \ n), 

Q 



n=0 



where 



CQ = Cn= Y I ^« I ^s) '^Q'. 



(20) 



(21) 



where Q = {qi,q2,--- Qng) , and Q' = (qi, ■ ■ ■ ' ' ' • That is Q and Q' are 
equal except for the qubit acted upon by the one-body operator Qg- 



A better way to state the above result is to consider Eq. (21) for the case that 



n has Qs = and thus n ^ uq and to write out the sum over q'^ to get 

a„ = (0 Ifi.l 0)C„„ + (0 1)C„,, (22) 
where we introduced the partner to hq namely ni. For the case that n has 



Qs = 1 and thus n — )■ rii Eq. (21), with expansion of the sum over q'^ yields 

Cn, = {l\ns\0)Cn, + {l\ns\l)Cn,. (23) 

or written as a matrix equation we have for each uq, rti partner pair 



no 



This is not an unexpected result 



(0 1 1 0) (0 1 1 1) W a,o 
1 I fi, 1 0) (1 I 1 1) i 1 c„. 



(24) 



Equation (24) above shows how a 2 x 2 one-qubit operator Vts acting on qubit 



is changes the state amplitude for each value of ng. Here, tlq denotes a decimal 
number for a computational basis state with qubit ig having the value zero 
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and rii denotes its partner decimal number for a computational basis state 
with qubit ig having the value one. They are related by 

nl = r^o + 2"«~^^ (25) 

At times, we shall call 2""'"*^ the "stride" of the is qubit; it is the step in 
n needed to get to a partner. There are 2"'9/2 values of tiq and hence 2"'/2 



pairs no, rii. Equation (24) is applied to each of these pairs. In QCWAVE that 
process is included in the command Opl[^ 

Note that we have replaced the full 2"« x 2"« one qubit operator by a se- 
ries of 2"9/2 sparse matrices. Thus we do not have to store the full 2"« x 2"« 
but simply provide a 2 x 2 matrix for repeated use. Each application of the 
2x2 matrix involves distinct amplitude partners and therefore the set of 
2x2 operations can occur simultaneously and hence in parallel. That paral- 
lel advantage was employed in our QCMPI fortran version, using the MPI |9] 
protocol for inter-processor communication. The 7.0 & 8.0 versions of Mathe- 
matica include only master-slave communication and therefore this advantage 
is not generally available. It is possible to use MPI with Mathematica [lOJ, but 
only at considerable cost. Another promising idea is to use the "CLOJURAT- 
ICA" [TT] package, but that entails an additional language. So the full MPI 
advantage will have to wait until MPI becomes available hopefully on future 
Mathematica versions. 

In the next section, this procedure is generalized to two- and three-qubit op- 
erators, using the same concepts. 



3.2 Two-qubit operators 



The case of a two-qubit operator is a generalization of the steps discussed for 
a one-qubit operator. Nevertheless, it is worthwhile to present those details, 
as a guide to those who plan to use and perhaps extend QCWAVE. 

We now consider a general two-qubit operator that we assume acts on qubits 
and ig^, each of which ranges over the full 1, ■ ■ ■ , possible qubits. General 
two-qubit operators can be constructed from tensor products of two Pauli 
operators; we denote the general two-qubit operator as V. Consider the action 
of such an operator on the multi-qubit state | : 



^ Opl yields result of a one-body operator Q, acting on qubit "is"' in state ipo', the 
result is the final state ^pf. Called as: ipf = Opl[0, is, ■0o]- 
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V|vl/)„^=EC^Q VI Q) (26) 
Q 

= E--- E ■■• E C'q |gi)---(V|g.ig.2)) ■■■|g„,). 

91=0 93 1,982=0 9na=0 



Here V is assumed to act only on the two g^i, qs2 qubits. The (V| g^i qs2)) term 
can be expressed as 

1 

V\qsiqs2)= E Usl 9s2)(?sl |V| gsi gs2) (27) 

9^1.9^2=0 



using the closure property of the two-qubit product states. Thus Eq. (27) 
becomes 



v|^k=EcQV|Q)= E--- E ■■■ E ■■■ E E (28) 

Q gi=0 9sl=0 9a2=0 q„g=0 q'^^,q'^2=0 

Cq {q'si(l's2\^\ Qsiqsi) \ qi)--- \ q'sM ' ' ' I 

Now we can interchange the labels qsi qsi,qs2 ^ qs2 use the label 
Q to obtain the algebraic result for the action of a two-qubit operator on a 
multi-qubit state, 



2"i-l 

V|^)n, = E^Q IQ)= E Cn \ n), (29) 

Q n=0 

where 

1 

CQ = Cn= E {(lsiqs2 I I q'sM Cq', (30) 
9^1-9^2=0 

where Q = (^qi, q2,--- qn,) , and Q' = (^qi, ■ ■ ■ g^i " " " " " " • That is Q and 
Q' are equal except for the qubits acted upon by the two-body operator V. 



A better way to state the above result is to consider Eq. (30) for the following 
four choices 



^00 -> {qi ■ 


■qsi 


= 0- 


■qs2 


= 0,- 


■■gnj 


noi (gi ■ 


■ qsi 


= 0- 


■ qs2 


= 1,- 


■■(In,) 


nw-^{qi ■ 


■qsi 


= 1- 


■qs2 


= 0,- 


■■In,) 


nil -> (gi • 


■qsi 


= 1- 


■qs2 


= 1,- 


■■In,) 



(31) 

where the computational basis state label n^^^^q^^ denotes the four decimal 
numbers corresponding to Q = {qir ■ ■ qsi ■ ■ ■ qs2 ■ ■ ■ qn„) ■ 
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Evaluating Eq. (30) for the four choices Eq. (31) and completing the sums 
over q'si,q'g2i the effect of a general two-qubit operator on a multi-qubit state 
amplitudes is given by a 4 x 4 matrix 



where V, 



ij-M 



(c \ 




VoO;00 


VoO;01 


VoO;10 


VoO;ll ^ 




(c \ 






Voi;00 


Voi;01 


Vol; 10 


Voi:ll 




c 


n 




VlO;00 


VlO;01 


VlO;10 


Vl0;ll 




r 

"-'"10 


\ ^"11 / 




\ Vll;00 


Vll;01 


Vll;10 


Vll;ll ) 






= {hi 


1 V 


1 Kl). 


Equat 


ion (32) shows how a 4 



(32) 



operator V acting on qubits i^i, '^s2 changes the state amplitude for each value 
of noo- Here, noo denotes a decimal number for a computational basis state 
with qubits is\i"i's2 both having the values zero and its three partner decimal 
numbers for a computational basis state with qubits isi,is2 having the val- 
ues (0, 1), (1,0) and (1, 1), respectively. The four partners noo, ''^-oi, ''^■lo, ""-ii, or 
"amplitude quartet" , coupled by the two-qubit operator are related by: 

noi = noo+2"'-'^^ nio = noo+2"'"'=' rin = noo+2"'-'=^+2"«-*=^ (33) 

where is2, is2 label the quarks that are acted on by the two-qubit operator. 

There are 2""'/4 values of rioo and hence 2"«/4 amplitude quartets noo, "^oi, ""-lo, '^ii- 



Equation (32) is applied to each of these quartets for a given pair of struck 
qubits. In QCWAVE that process is included in the command Op2p^ 

In this treatment, we are essentially replacing a large sparse matrix, by a set 
of 2"''/4 4x4 matrix actions, thereby saving the storage of that large matrix. 



3.3 Three-qubit operators 



The above procedure can be extended to the case of three-qubit operators. 
Instead of pairs or quartets of states that are modified, we now have an octet 
of states modified by the three-qubit operator and 2""'/8 repeats to cover the 
full change induced in the amplitude coefficients. For brevity we omit the 
derivation. In QCWAVE that process has been implemented by the command 

OpsQ 

^ Op2 yields result of a two-body operator U acting on qubits "is" and "is2" in 
state ipo; the result is the final state V'/- Called as: tpj = Op2['D,,isl,is2,ipQ]. 
® Op3 yields result of a three-body operator $7 acting on qubits "zsl","is2" 
and "is3" in state -00 ; the result is the final state ipf. Called as: ipf = 
Op3[0, isl,is2, zs3, V'o]- 
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4 IMPLEMENTATION OF GATES ON STATES 



We now present some sample cases in which the above gates are apphed to 
state vectors. 



4-1 One-qubit operators 

Consider a state vector for Ug = 3 qubits defined by |\l/)3 =| 000), which in 
vector form is 



fr \ 




M 

























— > 

























[Cj) 







Now have a Hadamard act on qubit 1 by use of the command Opl[H, 1, "^s]- 
The result is displayed in vector form and then in Dirac form by use of the 
DForm command in Figure [Tj 

One can act with Hadamards on every qubit, by either repeated use of Opl 
Figure |2| or by the command f2ALL['H,psi], which is illustrated in Figure |3} A 
Dirac type notation is also available as illustrated in Figure |4] 

In QCWave.m the command Opl is given as a Module see Figure |5| which 
makes use of the command Pickl, Pickl selects the pairs of decimal labels 
which differ only in the "is" qubits value of 1 and 0. Then all such pairs are 
swept through. Examples of Pickl are presented in the Tutorial. 



4-2 Two- qubit operators 

The typical two-qubit operators are the CNOT, and controlled phase oper- 
ators. General two-qubit operators can be constructed from tensor products 
of two Pauli operators, as discussed earlier. An example from QCWave of 
application of a CNOT gate is presented in Figure [6j 
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Out[98]= 



ln[98]:= psi 

1 ' 







I } 



ln[99]:= Opl ["K, 1, psi] 



1 






1 







Out[99]= 



ln[100]:= DForm[Ss] 

1 

Out[100]= + I 000 > + 



I ICQ > 



Fig. 1. One Hadamard example. Here psi=| 000). 
A Dirac type notation is also available as illustrated in Figure [7] 

In QCWave.m the command Op2 is given as a Module see Figure |8| which 
makes use of the command Pick2, Pick2 selects the quartet of decimal labels 
which differ only in the "isl" and "is2" qubit's values of 1 and 0. Then all such 
quartets are swept through. Examples of Pick2 are presented in the Tutorial. 



J^..3 Three-qubit operators 



The Op3 command is also provided in QCWave.m as a Module, which makes 
use of the command PickS, PickS selects the octet of decimal labels which 
differ only in the "isl," "is2," and "is3" qubit's values of 1 and 0. Then all 
such octets are swept through. Application to the TofoUi gate is provided in 
the Tutorial. 
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In[107]:= OplE'K, 3, Opl [-H, 2, Opl [-H, 1, psi] ] ] 



Out[107]= 



1 \ 

2 ^/2~ 
1 

2 4^ 
1 

2 

2 •TY 
1 

2 

1 

2 

1 

2 V2~ 
1 



2 V2~ / 

ln[108]:= DFonn[%] 



1 



Out[108]= 



2 V2" 



Oil > + 



I 000 > 



2 V2~ 



I 001 > + 



1 1 
I 010 > + 

2 V2~ 



2 V2~ 



I 100 > 



2 a/2~ 



I 101 > + 



2 V2~ 



2 V2 



I 110 > 



2 V2~ 



I 111 > 



Fig. 2. Hadamards on all three qubits example. Here psi=| 000). 



In[111]:= £JALL[W, psi] 

/ 1 

2 

1 

2 ^2~ 
1 

2 V2~ 
1 

2 

1 

2 

1 

2 V2~ 
1 

2 Vl" 

1 

2 V2~ 1 

ln[112]:= DForm[%] 



Out[111]= 



Out[112]= 



1111 
1000 > + 1001 > + 1010 > + I 

2^2" l-xpZ 2-\IY 2-{Y 



1 



Oil > + 



2 V2~ 



I 100 > + 



2 V2~ 



I 101 > 



2 a/2~ 



I 110 > + 



2 a/2~ 



I 111 > 



Fig. 3. Hadamards on all three qubits using the OALL command. Here psi=| 000). 
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In[113]:= psi 



Out[113]= 



1 







I I 



ln[114]:= XI = "K (1) I psi > ; 
X2 = "H (2) I XI >; 
X3 = -W (3) I X2 >; 

ln[117]:= DForm[Xl ] 
DForm [X2 ] 
DForm [X3 ] 

1 1 
Out[117]= + 1000 > + I 100 > 

V2~ V2~ 



Out[118]= 



Out[119]= + 



1111 

- 1000 > + - 1010 > + - |100 > + - |110 > 

2 2 2 2 



1 



2 V2" 
Oil > + 



1 1 
1000 > + 1001 > + 

2 ^/Y 2 V2~ 



I 010 > + 



2 V2~ 



I 100 > 



I 101 > 



2 V2~ 2 V2~ 2 V2~ 

ln[120]:= OplC-W, 3, Opl 2, Opl 1, psi]]] == X3 



I 110 > 



2 V2~ 



I 111 > 



Out[120]= True 



Fig. 4. Application of Hadamards in Dirac form. Here psi=| 000). 







{* One-body operator Q acts on qubit "is" in state "state"; 
Opl [ is, Q, state] ;Using QCMPI method*) 

Opl[ Q_?MatrixQ, is_, state_?VectorQ] :^ 
Module [ { wf , temp, il, jl, icas, nqu} , nqx-Log [2, Length [state] ] ; nqu-nqx; 
wf =state; temp=wf ; Do [ il=Pickl [is, icas, 1] ; jl=Pickl [is, icas, 2] ; 
temp[[il]]= n[[l]] [[l]].wf[[il]]+n[[l]] [[2]].wf[[jl]]; 
temp[[jl]]= n[[2]] [[l]].wf[[il]]+Q[[2]] [[2]].wf[[jl]], 
{icas, 1, 2^nqu/2 } ] ; temp] ; 





Fig. 5. The Opl command as stipulated in QCWave.m 
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In[132]:= Clear [chi] 
nqx = 3; 

chi = Flatten [ Ket [1] ® (Ket [0] (giKet [0] ) ] ; 

DForm [ chi ] 

Out[135]= + 1 I 100 > 
ln[136]:= CN = CN0T[2, 1,2] 

10 

10 

1 

I 1 j 



Out[136]= 



Apply CNOT to qubits 1 & 2 of state chi. 
View in Dirac Forms 

ln[137]:= Op2 [CN, 1, 2, chi] ; 
DForm [%] 

Out[138]= + 1 I 110 > 

Apply CNOT to qubits 1 & 3 of state ciii. 

In[139]:= Op2[CN, 1, 3, chi] ; 
DForm[%] 

Out[140]= + 1 I 101 > 

Fig. 6. Use of Op2 to apply CNOT gates. 





Apply CNOT to qubits 1 & 2, then 1& 3 of state chi. 


In[153]:= 


Yl = CN (1, 2) 1 chi >; 




Y2 = CN (1, 3) 1 Yl >; 


ln[155]:= 


DForm [ chi ] 




DForm[Yl] 




DForm[Y2] 


Out[155]= 


+1 1 100 > 


Out[156]= 


+1 |110 > 


Out[157]= 


+ 1 |111 > 



Fig. 7. Application of CNOT gates in Dirac form. 
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( * Two-body operator Q act s on qubit 

" is 1 " and " is2 " in state " st az e " ; 
Op2 [ isl, is2,Q, state] ;Using QCMPI method*) 

Op2 [Q_?MatrixQ, isl_, is2_, state_?VectorQ] := 

Module [ { wf , temp, i2,j2,k2, m2 , icase, nqu} , wf= state; teitip^wf ; 
nqx=Log [ 2 , Length [state] ] ; nqu=nqx; 
Do [ 

i2=Pick2 [ isl, is2, icase, 1] ; 
j2=Pick;2 [isl, is2, icase, 2] ; 
k;2=Pick2 [isl,is2, icase, 3] ; 
m2=Pick2 [isl, is2, icase, 4] ; 
tempi [12]]= a[ [1] ][[!]] .w£[ [12]] 

+n[ [1]] [ [2] ].wf [ [ j2] ]+Q[ [1] ] [ [3] ].w£[ [Ic2]]+B[ [1] ] [[4]] .wf [ [m2] ]; 
temp[ [ j2]]= n[ [2] ][[!]] ™£[ [12]] 

+n[[2]] [ [2] ].wf [ [ j2] ]+n[ [2] ] [ [3] ].w£[ [lc2]]+a[[2] ] [[4]] .w£[[m2] ]; 
temp[[k2]]= Q[ [3] ] [ [1] ] .»£[ [12]] 

+D[[3]] [[2]].wf [[j2]]+Q[[3]] [[3]].w£[[k2]]+Q[[3]] [[4]].wf[[m2]]i 
temp[[m2]]= Q [ [ 4 ] ] [ [ 1] ] .w£ [ [ 12 ] ] 

+0[[4]] [[2]].wf[[j2]]+Q[[4]][[3]].wf[[k2]]+Q[[4]] [ [ 4 ] ] .wf [ [m2 ] ] , 
{ lease, 1, 2 "nqu/ 4} ] ; temp] ; 



Fig. 8. The Op2 command as stipulated in QCWave.m. 
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5 THE MULTIVERSE APPROACH 



5.1 General remarks 



Mathematica 7.0 & 8.0 provide a master-slave parallel processing facility. This 
is not a full implementation of a parallel processing setup that allows commu- 
nication between the "slave" processors", such as used by the MPI |i9j protocol. 
If MPI were readily available in Mathematica then one could invoke the full 
capabilities discussed in QCMPI. That capability allows for the state vector 
to be distributed over several processors which increases the number of qubits 
that could be simulated. It is indeed possible to have Mathematica upgraded 
to include MPI slave to slave communication; as is available in the "POOCH" 
[To] package. However, since that is an expensive route and most Mathematica 
users do not have access to MPI, although one can hope for such a capability 
in the future, we have not invoked the full state distribution aspect. 

Nevertheless, the master-slave Mathematica 7.0 capability does provide for 
concurrent versions of a QCWave based algorithm to be run with different 
random error scenarios. Then an ensemble averaged density matrix can be 
formed which describes a real error prone QC setup. That opens the possibility 
of examining the role of errors and the efficacy of error correction methods 
using Mathematica. 

Therefore, we provide a sample of a parallel setup using some simple basic 
algorithms, where the parallel setup is described and explained in detail. The 
following steps are needed: (1) set up your Mathematica code to access several 
processors, see Appendix 1 for some help; (2) identify the processor number; 
(3) introduce random errors depending on the processor number; (4) assign 
a probability distribution for the various processors; (5) form an ensemble 
average over the processors and store that information as a density matrix on 
the master processor; (6) repeat these steps including, the algorithm, noise and 
finally error correction (EC) steps on all processors; (7) examine the resultant 
density matrix and its evolution to test the EC efficacy. This is an important 
program that we start by providing simple examples. 

A quantum system can evolve in many ways. Different dynamical evolutions 
are called paths [12] or histories [13]. We refer to these alternate evolutions 
as separate "universes" and a collection of such possibilities as a multiverse 
or ensemble of paths. Parallel processing provides a convenient method for 
describing such alternate paths. 
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5.2 The ideal and the noisy channels 



In our application, we assume that the main path follows an ideal algorithm 
exactly and the alternate paths incorporate the algorithm with possible noise. 
That noise is described by random one-qubit operators acting once, or with 
less likelihood twice. To describe this idea, which is realized in the notebooks 
MVl-Noise.nb, MV2-Noise.nb and MVn-Noise.nb, consider an initial 
density matrix po- For a pure state, po =\ ipo){i^o \, but it can be a general 
initial density matrix subject only to the conditions = p and Tr[p] = 1. 
The density matrix has 2^"' — 1 parameters and 2"' real eigenvalues A„ < 1 
with J2n = 1- The simplest one-qubit case has the form p = |(1 + P ■ a), 
where the real polarization vector P = Tr[(Tp], is within the "Bloch sphere". 



5.2.1 Storage case 

Consider a simple case where the ideal algorithm is simply leaving the state, 
as described by po? alone. This is a memory storage case. Ideally, p remains 
fixed in time. Assume however that this ideal case occurs with a probability 
p < 1 and that alternate evolutions occur with a probability e = 1 — p. For 
example, we take p = .8 and e = .2, corresponding to a 80% perfect storage 
and 20% possibility of noise. We also assume for more than 1 qubit cases that 
95% of the 20% noise (.2 x .95 — > 19%) involves a single one-qubit hit, while 
5% of the 20% noise (.2 x .05 — )■ 1%) involves two one-qubit hits. 

The ensemble average over all paths then yields a density matrix 

p^ = p Po + Q(i)t + Q(2)t . . . ^K) p^ ^K)t)^ (35) 

Hp 

where the operators VL^^\VL^'^'^ ■ ■ ■ Vl'^'^^'> act in each of the Up paths with a prob- 
ability Each of these Hp terms is evaluated on a separate processor, so that 
Up equals the total number of processors invoked. The above ensemble average 
preserves the trace: 

n=np 

Tr[py] = p Tr[po] + - E ^^[^^"^ Po ^^"^^^ = p + -n^ = 1. (36) 

n=l 

Here we assume that each fi(")tf](«) = and hence that Tr[l]('^) p^ fi(")t] = 

^ The two qubit case is of the form p = \{l + Pi-{a®l)+P2-{l®a) + ^ -(W ®'^)), 
where Pi = Tr[a ® 1)] & P2 = Tr[\ ® a] are the polarization vectors for qubits 1 
and 2 and C = Tr[^ ® is the 3x3 spin correlation tensor. Note the number of 
polarization plus correlations are 2^"'' — 1 = 3 (for one-qubit) and 15 (for two-qubits). 
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to 
n=l 

Po 

n=2 
n=np 

Fig. 9. The multiverse approach. Each curve represents a possible evohition of a 
quantum system of Ug qubits. The sohd curve represents the main evolution path, 
with probability p. The dotted curves represent the n*^ evolution path, with prob- 
ability The initial density matrix is po and the final pp- The unitary noise 

operators act on each path i. Here Q^'^ denote subsequent noise operators. The 
operators and denote the Ug algorithm operators which act on all paths. 

Tr[po ] = I. In addition, pj- = pf. 

This multiuniverse approach is illustrated in Figure [9j For the pure storage 
case the algorithm operators Qa, ■ ■ ■ are all set equal to unit operators. 
See later for a simple non-trivial case. 




5.2.2 Multiverse and POVM 

The above representation can also be cast in the POVM (Positive Operator 
Valued Measure) and Kraus operator form. The evolution can be expressed as 

n=np 

Pf=Y. Po (37) 

n=0 

where we define = ^ 1, and il"^'^ = y^fi". This is the form known 

as POVM, which can be deduced [H] from embedding a quantum system in 
an environment, which is then projected out. This evolution form can also be 
used to deduce the Lindblad [15] equation for the evolution of a density matrix 
subject to environmental interactions. Here we arrive at these forms from a 
simple multiuniverse approach. 



5.2.3 Multiverse and classical limit 

Our task is to set up this multiuniverse approach using the parallel, multi- 
processor features of Mathematica. Equation 35 describes the evolution of a 
density matrix after one set of operators act in the various possible paths. A 
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subsequent set of operators is described by 



n=n 

e 



PF = ppf + - E (^^"^ P/ ^^"^^)- (38) 



n=l 



As this evolution process continues to be subject to additional noise opera- 
tors, the density matrix evolves into a diagonal or classical form. In this way 
the noise yields a final classical density matrix with zero off-diagonal terms; 
this is the decoherence caused by a quantum system interacting with an en- 
vironment. If the qubit states are not degenerate and the noise is of thermal 
distribution, the density matrix in the classical limit will evolve towards the 
thermodynamic form exp(— ■^). At every stage, one can track the von Neu- 
mann entropy( S{p) = —Tr[p\n[p\\), the Purity ( Tr[p- p]), and the Fidelity ( 

^ of the system. In addition, the eigenvalues 
where in the classical limit the eigenvalues 
all approach and the entropy goes to S[p] — t- Ug. Subsystem entropy and 
eigenvalues can also be examined. 



F[p, Po] = Tr[^ ^ -p- ^]) 
of the system can be monitored 



5.3 Multiverse algorithms and errors 



The evolution of the density matrix, with noise included via the multiverse 
approach on the available processors, can also be implemented when an algo- 
rithm is included. The procedure is to act with the algorithm gate operators 
after each ensemble averaged density matrix is formed. The explicit expression 



is given in Equation 35 which for the nth step is 



Pn+l — 



Pn ^ Yl Ps^ksPn^is 



V k=l s=l,2 



(39) 



where Qj^ are the gates for the specific algorithm and flks are the noise op- 
erators on the kth processor for two cases, s = 1 denotes a one qubit noise 
operator hitting one qubit and s = 2 denotes one qubit operators hitting two 
separate qubits. The factors pk are assumed to be pi = .95 and p2 = .05, so 
that the one-qubit hits have a net probability (e * Pfc) of 19% and the double 
hit case a net probability of 1%. The algorithm operators Qa are applied to 
all all processors, but implemented by evaluation on the master processor. 

In QCWAVE, the above steps are implemented in the notebooks MVl-Noise, 
MV2-Noise and MVn-Noise, for systems consisting of 1,2, or n qubits. The 



key step is shown in Figure 10, where denE[n] denotes the ensemble averaged 



^ To evaluate this complicated expression, we find the eigenvalues of p.po and then 
form the sum a / | Aj | , to obtain a good approximate value. 
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density matrix at stage n, and the parallel part of the command distributes 
the evaluation of the noise over the "nprocs" processors, which is doubled to 
account for the "s" label in Equation 39 



denE[0] = p [initstate] ; 
denE[nv_] : denE[nv] = 
(OPA[nv] ) . (prob * denE [nv - 1] + eps * WaitAll [Siim [ParallelSubmit [ {i} , probx[i] * 

(OP[i, nv] .denE[nv- 1]) . (Adj[OP[i, nv] ] ) ], {i, 1, 2 * nprocs} ]]). (Adj [OPA [nv] ] ) 



Fig. 10. Here denE[0[] is the initial density matrix pQ and denE[nv] denotes the en- 
semble averaged density matrix after nv steps. The algorithm operators "OPA[nv]" 
are setup within the code, and the noise operators 0P[ i , nv ] act on the i^^ proces- 
sor at the nv^^ step. These noise operators are held fixed after they are randomly 
generated. 



A simple algorithm is illustrated in MVl-Noise; namely, one starts with the 
state I 0) which is then hit by a Hadamard H, and after an interlude of noise, 
another Hadamard hits, followed by a long sequence of noise. Without noise 
this process correspond to a rotation to the x-axis and then a rotation back 
to the z-axis. One also sees the polarization vector rotated, rotated back and 
then, after a sequence of noise hits, decay to zero and density matrix then 
evolves into a diagonal form, with both eigenvalues equal to 1/2. How is this 
simple process affected by noise during these steps? To answer that question 
the entropy, purity and fidelity evolution are tracked. The results from MVl- 



Noise are illustrated in Figure 11 



Another simple algorithm is illustrated in MV2-Noise; namely, one starts 

with the state | 00), qubit one is then hit by a Hadamard "H, and after an 

interlude of noise, a CNOT gate acts on both qubits. This is the algorithm for 

producing a Bell state CN0Ti,2?^i | 00) = ^(| 00)+ | 11)). This is followed 

by an inverse Bell operator T^iCNOTi 2, and then a long sequence of noise. 

The density matrix then evolves into a diagonal form, with all 4 eigenvalues 

equal to 1/4. The two polarizations, and the spin correlations are displayed 

along with the evolution of the eigenvalues, the entropy, purity and fidelity. 

7? 



The results from MV2-Noise are illustrated in Figures 12 



Clearly, more sophisticated algorithms can be invoked. We next consider how 
to monitor and correct for the noise. 



5.4 Multiverse and error correction 



5.4-1 Simulation of error correction 

Error correction (EC) typically involves encoding the qubits using extra qubits, 
then entangling those encoded qubits with auxiliary qubits. Measurements are 
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Fig. 11. The fidelity (left column), purity (central column) and entropy (right col- 
umn) evolution is displayed for one qubit initially in a state | Vo) =| 0), that is 
subject to a Hadamard at step 1 and another Hadamard at step 3, for various noise 
scenarios. The top plot is for X,Y and Z noise, the next plot down has Y and Z 
noise, and the next Z-noise only. The bottom plot is the noiseless channel, for which 
the fidelity returns to one, after dropping to l/\/2, which clearly reflects the action 
of 2 sequential Hadamards. 

made on the auxiliary qubits, so as not to disturb the original encoded qubits. 
Those measurements provide information as to whether an error has occurred, 
its nature and where it acted. Hence a remedial gate can be applied to undo 
the error. If desired, the encoded qubits can then be decoded and the origi- 
nal error-free qubit restored. That process is illustrated for simple X and Y 
errors on one qubit EC and for Shor's 9 qubit EC code in notebooks EC3x, 
EC3z and Shor9Tutorial. More sophisticated EC codes are available in the 
literature, along with a general theoretical framework [16]. This kind of EC 
has to be constantly invoked as an algorithm evolves, which is a rather awk- 
ward and qubit-costly process. Error in the gates themselves is an additional 
concern, usually one assumes perfect gates, with errors (noise) occurring only 
in-between application of the gates. 

For our purpose, instead of applying the procedures outlined in the above EC 
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Fig. 12. The fidelity (left column), purity (central column) and entropy (right col- 
umn) evolution is displayed for two qubits initially in a state | -00) =| 00), that is 
subject to a CN0Ti2'Hi and later a ^iCN0Ti2, for various noise scenarios. The 
top plot is for X,Y and Z noise, the next plot down has Y and Z noise, and the next 
Z-noise only. The bottom plot is the noiseless channel, for which the fidelity returns 
to one, after dropping to 1/ -v/2, which clearly reflects the action of the 2 sequential 
operators. 

notebooks, we simulate EC by a rather simple procedure. In the notebooks, 
the operators ss[i] are set equal to the Pauli operators s[i] to generate noise. 
By replacing ss[l] by s[0] (the 2x2 unit operator), all X- noise is turned off 
"by hand." Then a rerun is generated which has the same structure as with 
the noise, except the X-noise has been removed. Similar steps can be used to 
remove the Y and the Z-noise operators. In that way a set of results can be 
generated ranging from a full noise, to partial noise to no noise cases. Examples 



from MVl-noise and MV2-noise are presented in Figures 11 md 12 



If the user wishes to invoke other noise operators, that can be accommodated 
as well. For example, a general unitary random rotation can be used as a noise 



operator by invoking the form shown in Figure 13 Thus "ss[4] =UE1" is used 
to turn on such rotations. Replacing s[4] to s[0] again provides a way to turn 
this operator off to remove that noise element. 
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(* Wigner Rotation Function--used to represent a rotation with Euler angles a,b,c *) 

/ Cos [b / 2] -Sin[b / 2] \ 
UEl [a_, b_, c_] : = MatrixExpF- I a s [3] / 2] . .MatrixExp [ - I c s [3] / 2] ; 

\Sin[b/2] Cos[b/2] I 

(* Set ss[4] as the Unitary rotation matrix *) 

i4select = 1 

If [i4select - 0, 

ss[4] = UEl[RandoinReal[{0, 2 Pi}], RandomReal [ { 0, Pi}], RandomReal [ { , 2 Pi}]], 
ss[4] UEl [RandomReal [{0, 2 Pi}], RandomReal [{ , Pi}], RandomReal [{ , 2 Pi}]]] 

ss [4] 
ss [4] 

ss [4] == ss [4] 

If[i4select == 0, Print[" Fixed random unitary rotation matrix case"]. 
Print [ " Variable random unitary rotation matrix case" ] ] 



Fig. 13. This is an option to use a general rotation as the noise operator as an 
alternative to the Pauli operators. The s[4] location can be used to include this 
option. 

With this simple scheme, one can study many more noise and EC scenar- 
ios. For example, in the notebook MVn-Noise the case of an algorithm for 
5 qubits is presented. The algorithm consists of a Hadamard followed by a 
CNOT chain, i.e. CNOT15 CNOT14 CNOT13 CNOT12H1 | 00)., including 
noise in-between and after the 9th step. Detailed examination of the entropy, 
fidelity, purity and eigenvalues and noise is presented within that notebook. 
Of particular interest is the EC simulation results when the noise operators 
are turned of sequentially. 



6 ADDITIONAL FEATURES 



6.1 Amplitude displays 

The amplitude coefficients C„ can be displayed in various ways using the 
commands Amplitudes and MeterGraph as illustrated in Figure ?? 



6.2 Dime form 

The command DForm has already been demonstrated in Figures 1-4, and 6. 
Another Dirac form has be invoked in QCWAVE as shown in Figure 4. A 
more extensive Dirac notation scheme has been provided by Jose Luis Gomez- 
Muiioz et al. in Ref. [3j. 
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Fig. 14. Amplitudes command displays amplitudes as magnitude and phase bar 
graphs. -see QCWave.m 
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Fig. 15. MeterGraph displays amplitudes as argand plots-see QCWave.m 
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nqset = 2; 

initialKets [matr_] := initialKetsA [matr] 
Circuit [matTel2 , { } ] ; 

Show[%, PlotLabel -» Style [ " Two Qubit Teleportation ", 18], 

ImageSize -» 450] 
initialKets [matr_] := initialKetsq [matr] 



Two Qubit Teleportation 



l^2> 



|03> - H 



|04>4^ 

lOs) 

|06> 



H 



H 



® 2±r 



s — 

o 



l^i> 
\b2) 

l^3> 
l^4> 

l^6> 



Fig. 16. A sample circuit diagram as produced in notebook Teleporta- 
tionW is shown here. The initialKetsA sets the initial array sequence as: 
I gi>,| 92),! Oa)-- - I Oe). 

6.3 Circuit diagrams 



Illustrations of circuit drawing are included throughout the notebooks, with 
the CircuitTutorial notebook providing an overview. The commands are all 



defined in Circuits. m. One example is given in Figure 16 



Upgraded applications 



Upgraded versions of Grover [17j, Teleportation 
are included in the present version. 



and Shor [12] algorithms 
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7 CONCLUSION & FUTURE APPLICATIONS 



This package will hopefully be instructive and useful for applications to error 
correction studies. Hopefully users will contribute to improvements and exten- 
sions and for that purpose we are developing an interacting web page. When 
MPI becomes available on Mathematica, there will be another opportunity to 
upgrade QCWAVE to a full research tool. 

Application to explicit quantum computing problems, such as study of non- 
degenerate states and the associated phase factors, errors in gates themselves 
(where the gates are produced by explicit pulses), and the direct application 
of EC schemes are among the possible future applications. Novel EC schemes, 
such as stabilizing pulses or EC stable spaces could be additional fruitful 
applications. 
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A Setting up Processors with Mathematica 



In order to use parallel processing with Mathematica, one needs to first gain 
access to several processors. There are other ways to do this, but, the com- 
mands we used are given in Figure [XTTj Note that the user need to be sure that 
the ssh (secure shell) access is working and accesses the Mathematica com- 
mand on the other machines (which could be Macs or PCs or a combination 
of them). 



Preferences 



Interface Evaljation Appearance System Parallel Internet Connectivity Advanced 



Master Kernel Name: Local 



m 



Cenefal Preferences 



Launch parallel kernels: O Manual ©When needed OAtstartjp 

Evaluation failure handling: ©Retry Q Abandon 

21 Try to relaunch failed kernels Enable parallel monitoring tools 



Parallel Kernel Configuration 

Total number of configured kernels: 6 



r 



' Local Kernels Remote Services Cluster Integration Remote Kernels 

Remote Kernels lets you launch kernels on other computers using remote login. 



Add Host 



Remove Host 



Duplicate Host 



Kernels ^able Hostname: 



Loraines air mac 



5 \^{r3 a 
[^£3 □ 



MacBook 

^ Use custom launch command 



ssh -K -f ftaba kin el 0.0.1 

. 2 .'AppI icataon s.-'Mat h emati ca.ap p/to n terts/MacQ^ 

MathKernel -matlillnk - 

linkmode Connea -linkprococolTCPIP - 

linkname ■ ■ 2 ' ' -subkernel -noinit 



(^Disable Remote Kernels^ 



f Reset to Defaults") 



Help... ) (^Parallel Kernel Status...) 



Fig. A.l. The preferences setup to access otfier processors- use "Evaluation/Parallel 
Kernel Configuration/Remote Kernels" to get to tliis page. You also need to set the 
Local Kernels entry. The user also needs to establish an ssh link to the other ma- 
chines and be sure it properly accesses the Mathematica installed on all processors. 
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In addition, one needs to setup the basic programs and requisite packages on 
all the processors used. For that the initializations shown in Figure A.2 are 
needed. 



Needs [ "QDENSITY^ Qdensity ^ " ] ; 
ParallelNeeds [ "QDENSITY^ Qdensity ^ " ] ; 
Needs [ "QDENSITY^ QCWave^ " ] ; 
ParallelNeeds [ "QDENSITY^ QCWave^ " ] ; 
Needs [ "QDENSITY^ Circuits^ " ] ; 



Fig. A.2. The commands needed to invoke the packages on parallel processors. 
QCWave.m and Circuits. m are add-ons to the original QDensity.m package. 
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